/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAConvolve
  \brief Convolve the theoretical momentum distribution of the core to give an
  experimental distribution. Input: various exp resolutions, and momdis_th
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022/06/17
  \date Last modified: 2022/06/19 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <unistd.h>
#include "config.h"
#ifdef CATCH_CNOK
#include "catch.h"
#endif
#include "TARandom.h"
#include "TAHist.h"
#include "TAConvolve.h"
#include "TAAME.h"
#include "TAException.h"
#include "TAMath.h"
#include "TAInterpolate.h"
#include "TAIntegrate.h"
#include "TAParaManager.h"
#include "TAYaml.h"

#define polint TAInterpolateD::PolyInter
#define simp TAIntegrate<>::Simpson

static const TAAME *AME = TAAME::Instance(); // atomic mass evaluation table
static TARandom r; // random number generator
static constexpr double u0 = TAMath::uMeV(); // atomic mass unit in MeV
static constexpr double PI = TAMath::Pi();

TAConvolve::TAConvolve(int ap, int zp, int zc, double ek, const string &expres)
 : fAp(ap), fZp(zp), fZc(zc), fEk(ek){
  fMp = AME->massNUC(fZp, fAp);
  fMc = AME->massNUC(fZc, fAp-1);

  // calculate parallel momentum of the core in lab frame -- pcz //
  fPp = TAMath::EkPerUToP(fEk)*fMp/u0; // projectile lab momentum
  fBG0 = TAMath::BetaGamma(TAMath::EkPerUToBeta(fEk));

  // read experiment resolutions from the yaml config file //
  if(access(expres.c_str(), F_OK))
    TAException::Error("TAConvolve", "ctor: file %s non-existent", expres.c_str());
  const TAYaml v(expres);
  DBG = v["DBG"].as<double>();
  DCOSREL = v["DCOSREL"].as<double>();
  DBG0 = v["DBG0"].as<double>();
  SDBG0 = v["SDBG0"].as<double>();
  PDBG0 = sqrt(SDBG0*SDBG0 - DBG0*DBG0);
  fdPp = fPp * (PDBG0/fBG0); // sigma(fPp)
} // end the ctor

/// \retval exp momdis \param pvzcm pz of valence nucleon in proj. c.m. frame
double TAConvolve::ExpMom(double pvzcm) const{
  // gaus convolution of the projectile momentum
  const double pp = r.Gaus(fPp, fdPp);
  // update projectile beta and gamma
  const double b0 = TAMath::PToBeta(pp, fMp), g0 = TAMath::Gamma(b0);
  // sample cos(theta) for pvzcm so as to obtain total pvcm
  // const double gccm =
    // TAMath::Gamma(TAMath::PToBeta(pvzcm/cos(r.Uniform(0., PI)),fMc));
  // the Lorentz transform from projectile c.m. frame to lab frame
  // const double pcz = fMc/fMp*pp*gccm - g0*pvzcm;
  const double pcz = fMc/fMp*pp - g0*pvzcm;
  // convolve pcz with exp. resolution to give its measured value //
  // theoretically pc is needed, which is unavailable yet. We use pcz =~ pc
  // as cospl~0.99999, which only makes difference on sigma(pcz), not <pcz>
  const double bg = TAMath::BetaGamma(TAMath::PToBeta(pcz, fMc));
  // Apr. 2, 2023: 2* -> the expMOMDIS is TA-NOTA, so the resolution doubles.x
  // Mar. 9, 2024: 2* removed -- the resolution does NOT change, but the counts do
  const double dpcz = pcz * sqrt(TAMath::Sum2(DBG/bg, DCOSREL));

  return r.Gaus(pcz, dpcz);
} // end member function ExpMom

/// generate random number with pdf as f(x) with length l
/// using hit or miss method
double TAConvolve::HitOrMiss(int l, const double *x, const double *f){
  const double a = x[0], b = x[l-1]; // [a,b]
  double c = 0., r1 = 0.;
  for(int i = l; i--;) if(c < f[i]) c = f[i]; // max[f(x)]

  while(true){
    r1 = r.Uniform(a, b);
    if(r.Uniform(c) <= polint(x, f, l, 3, r1)) break; // hit
  } // end while

  return r1;
} // end member function HitOrMiss

// cc: normalization constant
void TAConvolve::ExpMomDis(int l, const double *x, const double *f, double cc,
    vector<double> &xconv, vector<double> &fconv, ofstream &fout) const{
  const int n = 100;
  const double pc0 = fBG0*fMc;
  TAH1D h("h", "h", n, pc0-500., pc0+500.);
  // generate different random number sequence each time of calling
  r.SetSeed(rand());
  for(int i = 5000000; i--;){
    // if(i%10000 == 0) cout << fPp << endl;
    h.Fill(ExpMom(HitOrMiss(l,x,f))); // convolve
  } // end for over i
  h.Scale(cc/h.Integral("width")); // normalize to cc

  // output the convoluted momdis to array and file //
  xconv.clear(); fconv.clear();
  fout << std::fixed << std::setprecision(4);
  for(int i = 0; i < n; i++){
    xconv.push_back(h.GetBinCenter(i+1));
    fconv.push_back(h.GetBinContent(i+1));
    fout << xconv[i] << " " << fconv[i] << endl;
  } // end for over i
  fout.close();
} // end member function ExpMomDis

#ifdef CATCH_CNOK
TEST_CASE("Histogram and random numbers", "[hist]"){
  SECTION("general"){
    TAH1D h("h", "h", 3, -5, 25);
    h.Fill(-10);
    h.Fill(0); h.Fill(0);
    h.Fill(10); h.Fill(20);
    h.Fill(30); h.Fill(30); h.Fill(40);

    // h.Print();
    CHECK(h.GetBinContent(0) == 1);
    CHECK(h.GetBinContent(1) == 2);
    CHECK(h.GetBinContent(2) == 1);
    CHECK(h.GetBinContent(3) == 1);
    CHECK(h.GetBinContent(4) == 3);
    CHECK(h.GetBinCenter(2) == 10);
  } // end section
  SECTION("gaus"){
    TAH1D h("h", "h", 30, -4., 4.);
    for(int i = 10000; i--;) h.Fill(r.Gaus());
    h.Scale(0.05);
    h.Print();
    const int a = h.GetBinContent(15), b = h.GetBinContent(16);
    CHECK(double(fabs(a-b)/a) < 0.2);
    CHECK(h.GetBinContent(30) <= 1);
  } // end section gaus
  SECTION("ufm"){
    TAH1D h("h", "h", 10, 0., 10.);
    for(int i = 10000; i--;) h.Fill(r.Uniform(10.));
    h.Scale(0.05);
    h.Print();
    const int a = h.GetBinContent(1), b = h.GetBinContent(8);
    CHECK(double(fabs(a-b)/a) < 0.1);
  } // end section gaus
} // end TEST_CASE

// TEST_CASE("Random Number Generator", "[rdm]"){
//   TARandom r;
// } // end TEST_CASE
#endif
